// Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.

using System;
using dnAnalytics.Resources;

namespace dnAnalytics.LinearAlgebra.Solvers
{
    /// <summary>
    /// The convergence monitor that controls the iterative solvers.
    /// </summary>
    /// <remarks>
    /// The convergence check algorithm was taken from: <br/>
    /// Templates for the solution of linear systems: Building blocks
    /// for iterative methods
    /// <br/>
    /// Richard Barrett, Michael Berry, Tony F. Chan, James Demmel,
    /// June M. Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo,
    /// Charles Romine and Henk van der Vorst
    /// <br/>
    /// Url: <a href="http://www.netlib.org/templates/Templates.html">http://www.netlib.org/templates/Templates.html</a>
    /// <br/>
    /// Algorithm is described in Chapter 4, section 4.1, page 57
    /// </remarks>
    public sealed class ConvergenceMonitor
    {
        /// <summary>
        /// The infinity norm of the matrix on which the solver is running.
        /// </summary>
        private readonly double mMatrixNorm = double.NaN;
        /// <summary>
        /// The maximum (relative) value for the residual that must be achieved before convergence can
        /// be declared.
        /// </summary>
        private readonly double mMaximumResidual = 1e-10;
        /// <summary>
        /// The maximum number of iterations that may be performed before
        /// stopping the matrix solver.
        /// </summary>
        private readonly int mMaximumIterationCount = 100;

        private double mVectorNorm = double.NaN;

        private SolutionStatus mStatus = SolutionStatus.None;

        /// <summary>
        /// Initializes a new instance of the <b>ConvergenceMonitor</b> class.
        /// </summary>
        /// <param name="matrixNorm">The matrix norm on which the convergence calculations are based.</param>
        /// <remarks>The minimum residual is set to 1e-10 and the maximum number of iterations is set to 100.</remarks>
        public ConvergenceMonitor(double matrixNorm)
        {
            mMatrixNorm = matrixNorm;
        }

        /// <summary>
        /// Initializes a new instance of the <b>ConvergenceMonitor</b> class.
        /// </summary>
        /// <param name="maximumNumberOfIterations">The maximum number of iterations the <see cref="IIterativeMatrixSolver"/> is allowed to perform.</param>
        /// <param name="maximumResidual">
        /// The maximum (relative) value for the residual that must be achieved before 
        /// convergence can be declared.
        /// </param>
        /// <param name="matrixNorm">The matrix norm on which the convergence calculations are based.</param>
        /// <include file='../../../../examples/trunk/LinearAlgebra/Solvers/ConvergenceMonitorConstructor.xml' path='example'/> 
        public ConvergenceMonitor(int maximumNumberOfIterations, double maximumResidual, double matrixNorm)
        {
            mMaximumIterationCount = maximumNumberOfIterations;
            mMaximumResidual = maximumResidual;

            mMatrixNorm = matrixNorm;
        }

        /// <summary>
        /// Stores the vector norm of the right hand side vector.
        /// </summary>
        /// <remarks>
        /// The vector norm of the right hand side vector is used to calculate
        /// the stop criteria for the iteration sequence. This vector norm must be
        /// set before the actual iteration starts.
        /// </remarks>
        /// <param name="norm">The norm.</param>
        public void VectorNorm(double norm)
        {
            mVectorNorm = norm;
        }

        /// <summary>
        /// Indicates if the iterative process should continue or not. More information
        /// about the reason for stopping the process can be obtained by getting
        /// the <see cref="SolutionStatus"/> through the <see cref="ConvergenceMonitor.SolutionStatus"/>
        /// property.
        /// </summary>
        /// <param name="currentIterationNumber">The current iteration number.</param>
        /// <param name="solutionVector">The current solution vector.</param>
        /// <param name="residualVector">The current residual vector.</param>
        /// <returns><b>true</b> if the iterative process may continue; otherwise, <b>false</b>.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="solutionVector"/> is <b>null</b>.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="residualVector"/> is <b>null</b>.</exception>
        /// <include file='../../../../examples/trunk/LinearAlgebra/Solvers/ConvergenceMonitorShouldCalculationContinue.xml' path='example'/> 
        public bool ShouldCalculationContinue(int currentIterationNumber, Vector solutionVector, Vector residualVector)
        {
            if (solutionVector == null)
            {
                throw new ArgumentNullException("solutionVector", Strings.NullParameterException);
            }
            if (residualVector == null)
            {
                throw new ArgumentNullException("residualVector", Strings.NullParameterException);
            }

            // check the iteration count. If we are over the maximum number of
            // iterations we stop
            if (currentIterationNumber >= mMaximumIterationCount)
            {
                mStatus = SolutionStatus.IterationBoundsReached;
                return false;
            }

            // Store the infinity norms of both the solution and residual vectors
            // These values will be used to calculate the relative drop in residuals
            // later on.
            double residualNorm = residualVector.InfinityNorm();

            // Check the residuals by calculating:
            // ||r_i|| <= stop_tol * (||A|| * ||x_i|| + ||b||)
            double stopCriterium = ComputeStopCriterium(solutionVector.InfinityNorm());

            // First check that we have real numbers not NaN's.
            // NaN's can occur when the iterative process diverges so we
            // stop if that is the case.
            if ((Double.IsNaN(stopCriterium)) || (double.IsNaN(residualNorm)))
            {
                mStatus = SolutionStatus.Diverged;
                return false;
            }

            // Compute the relative epsilon
            double epsilon = mMaximumResidual * 1e-3;
            // Note the substraction can be negative but the stopCriterium will never be
            // So if the substraction is negative then the residualnorm is actually smaller
            // than the stopCriterium and we can safely return.
            // The reason for not putting it all in a System.Math.Abs is that if the residual is
            // zero (0) then the check would fail if we had a Abs in there.
            if ((System.Math.Abs(residualNorm) - stopCriterium) <= epsilon)
            {
                mStatus = SolutionStatus.Converged;
                return false;
            }

            // Set the status to none. We're still running so there is no
            // real status.
            mStatus = Solvers.SolutionStatus.None;
            return true;
        }

        private double ComputeStopCriterium(double solutionNorm)
        {
            // This is criterium 1 from Templates for the solution of linear systems.
            // The problem with this criterium is that it's not limiting enough. For now 
            // we won't use it. Later on we might get back to it.
            //return mMaximumResidual * (System.Math.Abs(mMatrixNorm) * System.Math.Abs(solutionNorm) + System.Math.Abs(mVectorNorm));
            
            // For now use criterium 2 from Templates for the solution of linear systems. See page 60.
            return mMaximumResidual * System.Math.Abs(mVectorNorm);
        }

        /// <summary>
        /// Returns the status of the solution. Note that this is only of
        /// use once the iterative process has been stopped.
        /// </summary>
        /// <value>The <see cref="SolutionStatus"/>.</value>
        /// <include file='../../../../examples/trunk/LinearAlgebra/Solvers/ConvergenceMonitorSolutionStatus.xml' path='example'/> 
        public SolutionStatus SolutionStatus
        {
            get
            {
                return mStatus;
            }
        }
    }
}
